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Abstract 

We study the relaxation of force distributions in the q-model, assuming a uni- 
form q-distribution. We show that "diffusion of correlations" makes this relax- 
ation very slow. On a d-dimensional lattice, the asymptotic state is approached 
as /( 1_d )/ 2 ! where I is the number of layers from the top. Furthermore, we de- 
rive asymptotic modes of decay, along which an arbitrary short-range correlated 
initial distribution will decay towards the stationary state. 

PACS numbers: 02.50.Ey, 81.05.Rm 



1 Introduction 

To answer basic questions with simple models is the favorite topic in the work of Bob 
Dorfman. This contribution, dedicated to him on the occasion of his 65th anniversary, 
concerns the force relaxation in a model for granular media. The model is certainly 
simple but whether the question we answer is a basic one, we leave to his judgment. 

The first idealization of granular media is to consider them as a pack of equally sized 
beads. The second idealization is to put the beads on a regular lattice, with periodic 
boundary conditions horizontally in order to avoid boundary effects. The problem is 
the propagation of the force exerted on the top layer downwards to the deep lower 
layers of the bead pack. Now it might seem that all randomness, characteristic for 
granular media, has been lost in this idealization. However the forces which the beads 
exert on each other are supposed to be transmitted in a random way. Let be the 
force in the downward direction on the ith bead in a layer. This bead makes contact 
with a number of z beads in the layer below, which we indicate by the indices i + a. 
The a's are displacement vectors in the lower layer as shown in Fig. |I] (a) for a 2- 
dimensional triangular packing and in (b) for a 3-dimensional fee packing. Bead i 
transmits a fraction qf of the force /, to the bead i + a underneath it. The random 
element in the model is that the fractions qf are taken stochastically from a uniform 
distribution satisfying the constraint 

£<tf = l- (1) 

a 

This is the q-model as introduced by Liu et al. The problem we consider is: given 
the force distribution P (fi, ■ ■ ■ , f N ) in the top layer, what is the force distribution 
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in the lower layers and in particular how does the distribution approach its asymp- 
totic value? Recent studies |2], U addressed the relaxation of the second moments for 
general q distributions on a triangular lattice. Restricting ourselves to the uniform q 
distribution, we investigate the relaxation of the full distribution for arbitrary lattices. 

In order to write the equation for the force distribution we consider the force /j on 
the jth bead in layer. It is composed of forces fj- a from the layer on top of it as 

fj = E fi-Jj-a- (2) 
a 

We have left out the so-called injection term representing the weight of the particles.^] 
As a consequence of the regularity of the lattice, a bead is supported by z beads in the 
next layer and reciprocally a bead in the lower layer also has z beads pressing on it. 
Note that, due to the fluctuations in the g's, the sum over the q^_ a in (H) does not add 
up to 1, which means that even a set of equal fj- a will lead to fluctuations in the /j. 




Figure 1: The displacement vectors a in the q-model for (a) the triangular packing 
(side view) and (b) the fee packing (top view) 

The transformation of the force distribution from one layer to the next layer can 
be written as 

P'(f') = dff DqH8(f^-^2 q f_J^ a )P(f), (3) 

j a 

— * 

where we have introduced a vector notation for the forces in one layer / = (/i, • • • , /at), 
and for the integrations we use the abbreviations 

/ d f = II jf dfi , jDq = HD qi = l[ (z-l)\ J[ jf* dq? 6 if ~ lj ■ (4) 

As one sees from the integration over the g", we have taken a uniform and normalized 
distribution over these stochastic variables. The restriction to a uniform g-distribution 
simplifies the mathematics while it still reproduces important experimental observa- 
tions [1. 

The recursive relation (|3[) between two successive layers contains the problem of 
force relaxation from top to bottom of the bead pack. For a number of properties it is 
more convenient to work with the Laplace transform of the distribution, 

P(s) = J dfeM-s- f) P(f)- (5) 

1 Leaving out the injection term is legitimate when the applied force is large compared to the 
particle weight or when the direction of propagation is perpendicular to gravity. 
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The Laplace transform of (|) leads to 

(6) 

a J 

The sum in the exponent can be rewritten as 

(7) 



P'(s) = J DqJ d/exp f- EE p (fi- 



E E 9?-a S -J fi-a = E E € S i+« ) fi 



by changing the summation variable j to i = j — a. So we can express the right hand 
side of (H) also in the Laplace transform and we get 



P'(s) = J DqP(s(q)), (8) 

with 

s i(o) = E^Si+a- (9) 

The projection of the full distribution P(f) to the single-site force distribution p(fi) 
is particularly simple in the Laplace language 

p( Si ) =P(0,---,0, Si ,0,---,0), (10) 

since = means a full integration over /j. 

Either (|^) or (^j) completely defines the propagation of the forces for boundary 
conditions (N, M) , where N is the number of sites in a layer and M the number of 
layers, (both assumed to be large). For the ultimate asymptotic behavior one has 
M ^> N, but the N ^> M case is physically more relevant and therefore the main 
focus of this paper. It is also not sensitive to the (periodic) boundary conditions that 
we have chosen. 

Coppersmith et al. p derived a stationary state of ©, 

P*(f) = IM/i) » P*(/) = z *f~ l exp(-zf)/(z - 1)!. (11) 

It is a product state without any correlation between sites. They provided numerical 
evidence that it is indeed the asymptotic force distribution. This paper concerns the 
decay towards the stationary state, which turns out to be algebraic rather than ex- 
ponential, as is implicit in ||. We will show how "diffusion of correlations" accounts 
for this slow relaxation. Furthermore, we will pay attention to the question of which 



initial distributions evolve to the stationary state (11) and in what sense the approach 
has to be understood. 

The first part of this paper deals with the decay of the distribution on the level of 
its moments. We start out by showing that the distribution of the total force in a layer 
is invariant under the recursion. Then we illustrate the propagation of forces in the 
system by considering the evolution of the first few moments of the force distribution. 
In the second part of this paper we construct solutions of the recursion relations, which 
dominate the asymptotic decay towards the stationary state (|TT|), with an emphasis 
on the relaxation of the single-site force distribution. In the concluding remarks we 
comment on the results that we have obtained and discuss which of these can be 
generalized to arbitrary q distributions. 
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2 The distribution of total force 



There is one obvious invariant in the problem: the total force on the particles of one 
layer 

F = Ef<- ( 12 ) 

i 

This can be seen by summing (E3) over all j 

by again changing the summation variable j to i = j — a. The magnitude of F is 
irrelevant for the problem since, due to the linearity of the force law @, an overall 
scaling of the forces is possible without changing the physics of the problem. Thus 
we could fix the value of F and only consider distributions having strictly this value, 
but that is mathematically not very convenient. The distribution of the total force is 
obtained from the distribution of the forces as 

R{F)= fdfS(F-J2fi)P(f)- (14) 

i 

Since relation fll3|) holds for any set of g's, one finds from (|j) that R(F) is invariant. 

So the initial distribution R^ (F) dictates that of the asymptotic state. This seems 
a strong restriction on the relaxation of the distribution function, but in practice it is 
rather a warning on what quantities to inspect. As an example consider the total force 
distribution of the stationary state (^lj), which reads 

R*(F) = Z N F zN ~ l exp(-zF)/(zN (15) 

It is a sharp distribution with a mean and variance 

(F) = N, (F 2 )-(F) 2 = N/z. (16) 

Note that fixing the total force F excludes ([H]) as the stationary state! The distri- 
bution ( jl5|) is very reminiscent of the energy distribution in the canonical ensemble. 
Indeed this analogy, also mentioned in ||, is illuminating. For large N the total force 
distribution (|15|) is narrow, in the same way as the total energy distribution in the 
canonical ensemble. Now it is well known that the correlation functions of the canoni- 
cal ensemble and the micro-canonical ensemble (in which the energy is fixed) , coincide 
for all distances (to order 1/iV). But integrals over the correlation functions over all 
distances differ in the two ensembles. A similar subtlety arises here. As long as we 
consider spatial correlations between forces, the distribution of the total force is unim- 
portant, if it is sharp in the sense of (0). However, for sums over the correlations, 
differences will appear. Consider for instance the force correlations between two sites, 
defined as 

(fif i+n ) = Jdffif i+n P(f). (17) 
Summing (fifi+ n ) over i and n yields 

J2( fifi+n) = J df X fifi+n P(f) = (F 2 }. (18) 
i,n ' i,n 

If F is allowed to fluctuate, this sum clearly differs from the result obtained for fixed 
values of F. The fluctuations which we allow in the total force are limited to 

(F 2 ) - (F) 2 = cN, (19) 

with c of order 1. Then they are unimportant in the thermodynamical limit. 
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3 Diffusion of the first moment 



Forces applied to one bead in the top layer will diffuse into a wider region in the 
lower layers. To see this, consider the first moments ( fj ) of the distribution for an 
inhomogeneous initial distribution. The relaxation of these moments is given by 

( Si )' = I dff Dq £ q«_ a f 3 _ a P(f) = f DqJ2 fj-a )• (20) 

Integrals over the g"s are elementary, but this one is trivial as each gj_ a gives the same 
answer, and with (M it must be equal to 1/z. So the result becomes 

</;>' = \ E</i-«)- ( 21 ) 

z a 

This simple recursion law tells us that a distribution with only a non-vanishing moment 
on one site in the top layer gradually spreads over the lower layers like a diffusion 
process. 

In order to understand the role of system size it is useful to inspect the fourier 
transform 

m(k) = X;</<)exp(ik.r i ), (22) 

i 

which relaxes from layer to layer as 

m'(k) = A(k)m(k), (23) 

with the transmission function 

A(k) = -£exp(ik-r a ). (24) 

z a 

In accordance with the conservation of the total force we see that the mode k = is 
conserved as A(0) = 1. All the other modes decay exponentially since |A(k ^ 0)| < 1. 
So we find asymptotically 

m (oo) (k) = 5 kiO m(0), (25) 
with m(0) the average of the (initial) total force. Translating it back to space yields 

(f i )^ = (F)/N } (26) 

i.e. all sites feel the same average force regardless of the initial distribution. 

Clearly this result holds in the limit M 3> N. In the other limit IV>Mwe may 
replace the summation over k by an integration over the Brillouin Zone (BZ) of the 
layer and one gets the standard diffusion process from the integral 

(fj) {M) = 77— [ dkm(k)exp[Mlog(A(k))-*k •!•,]. (27) 
Vbz j bz 

In the limit of large M the integral can be evaluated by the saddle point method, 
using the fact that the dominant contribution comes from small k. This allows us to 
approximate A(k) by 

log[ A(k) ] ~ log[ 1 - £(k • r„) 2 /2 ] ~ —D k 2 , (28) 
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which renders the integration effectively a Gaussian. The diffusion constant D depends 
on the latttice structure of the layer, with D = 1/8 for the triangular packing and 
D = 1/12 for the fee packing. Thus one sees that the ratio of M/N decides whether 
we should take a discrete sum over the k variables (M ^> N) or whether we should 
integrate over the Brillouin Zone (N 3> M). In the first case, only k = matters and 
we have diffusion crossing the periodic boundary conditions in the layer. In second 
case, a region around k = has an influence and there are no effects from the finite 
size of the layer. 

Having seen the washing out of spatial inhomogeneities by diffusion, we now con- 
centrate on translational invariant initial conditions and on the case N ^> M. 

4 Evolution of the second moment 

As we are free to put ( fi ) = 1, the relaxational properties show up only in the higher 
moments of the force distribution. We therefore examine how these moments evolve 
under the recursion. The transformation of second moments ( fifi+ n ) is obtained by 
combining (|T7|) and (EH as 

( fjfj+n >' = E (/ D ?$-a ^+n-S) ( fj-afj+n-a' )■ (29) 
a,a' 

The g-integrals can be worked out and yield 

Dq g£_ a qi+n-a' = 4 + ( --27TTT\ + \ ) S o,n+*-*' ■ (30) 



z 2 \ z 2 {z + l) z{z + l 
Inserting this into (|29[ ) gives the following recursion scheme 



fjfj+n )' — 72 £( fj-afj+n-a' )+9n{fj)- (31) 



a, a. 



The function g n incorporates all corrections due to the Kronecker 5's in (30) and the 
only non-zero terms are 

z-l 1 

9o — —, — ! 7T i 9-i — — T7 — ~~rr- (32) 

z\z + 1) z z {z + 1) 

The index 7 denotes a nearest neighbor position. It is easily checked that the difference 
of two unequal displacement vectors a and a' points to a nearest neigbor position. It 
is a matter of counting to verify that correlations of the type 

(/i/i + n>* = C(l + ^ , n ) (33) 



are fixed points of relation (^TJ) for arbitrary C. This conclusion was also derived by 
Lewandowska et al. [[| for the specific case of a 2-dimensional triangular lattice and 
C — 1, (corresponding to (|ll|))- For C ^ 1 one has long-ranged correlations since, for 
large n, the average of the product does not approach the product of the averages (set 
equal to 1). C is related to the scale of the forces as one sees from (ITj3). This relation 



implies for (33) 



(F )* = C (N + N/z). (34) 
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For a sharp distribution of the total force in the sense of (|IED and normalization of the 
average (F) = N, there is little room for C to differ from 1, allowing only deviations 
of order 1/N. So we take C = 1 corresponding to the stationary state flTT|) and come 
back to the issue in the Appendix. 

To study the relaxation in more detail we consider the difference with respect to 
values of the stationary distribution 



An ( fifi+n ) ( fifi+n ) 



(35) 



which will of course relax according to the same relation (plf ). As an example, the flow 
diagram for the A n is depicted in Fig. ^| for the triangular lattice. The advantage of 



-1 o 




t 



1/6 2/3 1/6 1/4 1/2 1/4 



Figure 2: The flow diagram for A$ in the triangular lattice, starting with only A °^ ^ 0. 
The horizontal index n indicates the relative distance. 

taking the difference is that its fourier transform is well behaved for distributions which 
have only short range correlations, since both the initial correlation function and the 
stationary values approach 1 for large distances n. The fourier transforms are defined 

as 



A(k) = VA n exp(zk- r n ), A n = — — [ dk A(k) exp(-zk • r n ), 

n v bz Jbz 

and A(k) transforms as 

A'(k) = \(k) \(-k) A(k) + A g(k). 
The correction function g(k) can be written as 

1 r 



g(k) = g + 9i exp(ik • r 



A(k)A(-k)]. 



(36) 



(37) 



(3* 



Note that A(0) is invariant, which implies 

n n 

as it should be according to relation (fT 



The recursion fl37|) can be solved by introducing the generating function 

i 



A(u,k) = Y, Ail) W 

1=0 



u 



(39) 



(40) 



where the sum extends over all layers I. Multiplying (|37|) by u l+1 and summing over / 
gives an algebraic equation for the generating function 



A(u, k) - A® (k) = u [ A(k) A(-k) A(u, k) + A (u) g(k) 



(41) 
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where the auxiliary function A (u) is defined as 



A (u) = J2< ) u l . (42) 
1=0 

A (u) is our prime interest as it contains the center coefficients Aq\ giving informa- 
tion about the relaxation of the single-site force distribution. The other A$ contain 



correlations between two sites. Using (|38|), we write the solution of (|4T|) as 



A(u, k) = 44tT7 ix + ~^ A M [ 1 - z I ■ (43) 

v ' 1 -«A(k) A(-k) z + 1 v ' \ l-uA(k)A(-k)y 

The first term will be shown to give purely diffusive behaviour. This can be understood 
from the recursion of ([H]): without the corrections g n , it becomes a regular random 
walk in a layer of dimension d — 1. The second term, which originates from the g n , 
contains a correction to the asymptotic amplitude and a term (proportional to 1 — u) 
representing a transient. Hence, we anticipate diffusive relaxation ~ /( 1_d )/ 2 . 
Using Eq. (|36p, the function A (u) satisfies 



A (u) = — dkA(u,k). (44) 
V B z Jbz 



Inserting expression ([4"3"D into (03) gives a closed equation for A (u) 



As we shall see, the large / behaviour is dominated by the behavior at small k. We 
have assumed the initial correlation to be short-ranged such that A^ (k) is a regular 
function of k at the origin. Therefore we may as well continue by considering first the 
uncorrelated initial distributions. Then only the value n = gives a contribution and 
we have initially 

A^(k) = A$\ (45) 

as depicted in Fig. ^|. A^ sets the overal amplitude of the relaxation process. 

Before discussing the asymptotics of the general case we begin by giving the results 
for the triangular packing with z = 2. In this case all the integrals can be carried out 
explicitly by writing them as contour integrals in the complex v = exp(ik) plane. For 
(HI) and (ID we get 

Ao(u) = + \a (u) (l - -±=Z£) . (46) 

This relation for Aq(u) can readily be solved as 

Aq{u) = ^ = I 3_\ = JA£^ f 1 _I v / r ^ + i(i_ M )_...V (47) 

V; s/I^l \2 + vT^/ 2VT^I V 2 V 4 V ; ) y ! 



We have made the expansion in powers of a/1 — u because each higher term leads to a 
weaker singularity and therefore to a faster decay. Tracing back these features to the 
form (^), we indeed see that the first term is responsible for the leading singularity, 
that the second modifies the amplitude by a factor 3/2 and that the third leads to higher 
powers and thus to transient behavior. The leading singularity gives an expansion 

1 ~ 1 (2m)! m . . 

y/l-u ( m! ) 2 
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Note that these are precisely the binomial coefficients one expects from the random 
walk described by the first term in fl3ll) . Using Stirlings formula for the factorials yields 
as the leading term for the A l 



A> ^ oo. ,49) 



(0) 



So we see that the singularity at u — 1 contains the asymptotic behaviour. This is 
also true for the general case. Defining the integral 

D{u) = -^f dk- 1 (50) 

Vbz Jbz 1 — u A(k)A(— k) 

we get for A (u) 

A (u) = 4° } D{u) + -^—A Q {u) [ 1 - (1 - u)D(u)}. (51) 

We again have not written the solution in closed form in order to see the origin of the 
terms. The first term contains the diffusive singularity at u — 1. This follows from 
the integral for D(u) which develops a singular denominator for u — 1 at the point 
k = 0. The type of singularity depends on the dimension d — 1 of the k integral and 
thus on the dimension of the system. The singularity can be obtained by the saddle 
point method using the approximation (|28|) for A(k) and reads (1 — up d ~ 3 >' 2 , where the 
power means a logarithm (for the fee packing). Consequently the decay is as /( 1_d )/ 2 
for large /. The other terms in ([51]) give a modification of the amplitude by a factor 
(z + l)/z and a faster decay term (due to a weaker singularity). 



In setting ( |45| ) we have made the restriction to initial states in which the second 
moments of the forces are uncorrelated. From the analysis it is clear that for the 
asymptotics it is a justified substitution for any initial state with short-ranged correla- 
tions. Thus the conclusion of this section is that for arbitrary initial distributions with 
short-ranged initial correlations the second moment of the force correlations relax to- 
wards the stationary distribution (jTTJ) . To be more precise, the approach is not perfect 
in view of the restriction (|39|). The values of all A^ saturate at a level ~ A^ /N when 
I starts to exceed the value N 2 ^ d ~ l \ Then the diffusion of the initial deviation goes 
around the periodic boundary conditions and the system will approach a state fl33|) 
with C — 1 of the order 1/N. This is the subtlety to which we alluded in Section |2|. 
The total sum over the A n corresponds to 

Y^A n = ((F 2 )-(F 2 y)/N } (52) 

n 

and when RS \F) is not perfectly equal to R*(F), the system will never fully relax to 
P*(f)- 

5 Stationary State and some Properties 

Coppersmith et al. || have determined a stationary state of Eq. (||), which is special 
in the sense that it is a product of single-site distributions. Since the proof of the 
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stationary distribution is vital for the decay of deviations, we briefly repeat its con- 
struction. It is based on the following identity which holds for uniform q distributions 
i 

D<li (l + AEatf&J = 5 1 + (53) 
Now if we apply this identity to the (normalised) distribution 

we see that each zth power of the factor splits up into z single powers of the neighboring 
sites. The product over all sites in a layer then restores the zth power of the denomi- 
nator, which makes P(s) an invariant. Putting the mean force equal to 1 requires that 
X = 1/z. So we have a stationary distribution 



Translating this back to the force distribution we get (|TTj) . 

Further information is deduced by differentiation of the identity ( |53|) with respect 
to A giving 



J q% (1 + XE a q?b a y+ l z { l } 1 + XbJ % 1 + Xb a , 



(56) 



This allows us to construct another invariant 

i 1 + ASi 

Inserting this expression into the recursion relation (|8|), we encounter the integral (|53"D 
for all sites, except for site i where we have to perform the integral of the left hand 



side of (|56|). That produces for a set of surrounding sites factors with the power z + 1 
and the numerics is such that adding them together yields again a sum as in (|57|) . It 
is not surprising that this additional invariant exists since we showed that ( |54| ) was 
invariant for arbitrary A or scale of the forces. Thus the invariance of ([57]) is nothing 
else than expressing this scale invariance of the recursion. In fact ( [5"TD follows directly 
from differentiating ([)4]) with respect to A. 

It is tempting and indeed rewarding to continue to differentiate the identity (p3|). 
Differentiation of (^) with respect to A yields 

D (E a q a b a ) 2 _ 1 / 1 \ Ml + <W)fr Q " r5g x 

(1 + A E a q a b a ) z+2 z{z + 1) 1 A, 1 1 + XbJ (1 + Xb al ) (1 + Xb a „) ■ 1 ' 



Now the distribution 



^(S) = -P*(S)E(-Z-) ' (59) 

; \Z + Si ' 



is not invariant since the right hand side of (58) produces terms with powers z + 1 of 
the denominator on two sites. We will show in next sections that (|59D leads to the 
"slowest mode" decaying towards to the stationary state. 

The fact that the distribution ([54] ) is stationary for arbitrary values of A leads to a 
multitude of other stationary distributions by differentiation of fl54|) with respect to A. 
The distribution (^) is an example. In the Appendix we comment on the relevance of 
these stationary states. From this section the main message is: 
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• In the transformation from layer to layer, the total power of the factors remains 
the same. This implies that we can select classes of perturbations which transform 
into themselves. 

• The sum of the coefficients of the terms is invariant under the transformation. 

6 The slowest mode 

Let us now investigate how a deviation from P*(s) will decay under recursion. We take 
an initial distribution of the type 

P(s) = P*{s)[l + AP{s)]. (60) 

In the spirit of the previous section we consider deviations of the form 

AP(s) = J2A n Q n (s), (61) 

n 

where the Q n {s) are products of two factors. 

'', ■ (62) 

It is clear from the previous Section that a AP(s) of the form (pi] ) is after the 
transformation again a sum over the Q n {s) 

AP\s) = J2A' n Q n (s). (63) 

n 

The transformed coefficients A' n are expressed in terms of the A n of the previous layer, 
using the formulae of the previous section, as 

<=^E A n+a-a> + 9n Aq. (64) 
^ a, a' 

On purpose we have chosen the same notation for the coefficients in the representation 
Q5"T| ) and in Section [| for the correlations in the forces. As one sees from comparing 
d&|) with d35|) and ([H]) the coefficients in both cases obey exactly the same recursion 
relation. This means that we can take over the conclusions of Section |j. The first one 
is that we can construct a stationary state of the form 

A* n = B{l + U nfl ). (65) 

The nature of this stationary state, having long-ranged correlations, will be discussed 
in the Appendix. The other relevant conclusion concerns the case where only Aq 7^ 
0. This is a minimal disturbance of the stationary state with only a modification 
of the single-site probability distribution and no correlations between sites.0 This is 
precisely the same initial condition as for the case of the second moment starting with 
an uncorrelated initial state. Thus the development of this deviation (|6T|) is exactly 
the same as that of the second moment of an arbitrary deviation from the stationary 
state. The difference is that we now have the development of the full force distribution. 

2 It is not a product distribution, but it can be seen as the linearization of a product distribution. 
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Projecting this joint distribution onto a single-site distribution, by putting all Sj = 
except for one, one finds 

^{8)=jr{s)(l+J$j^pJ. (66) 

So on the level of the single-site distribution the shape of the deviation remains in- 
variant, but the amplitude decays as Aq ~ £( 1-d )/ 2 . This means that not only the 
second moment, but the whole single-site distribution relaxes towards the stationary 
state (|TTD with a characteristic shape of the single-site distribution given by (|66|) . 



7 Higher Modes 

In this section we show that higher moments lead to slower modes. Let us start by 
looking at the third moment of the forces. The linearity of the force relation (0) 
guarantees that the third moments transform as a closed set 

( fjfj+nfj+m )' = ( / Dq qf_ a qf +n _ a , g" +m _ Q » J ( fj-a fj+n-a' fj+m-a" ) ■ (67) 

a,a',a" V J 

The q integral is complicated but elementary, with several exceptional cases due the 
equality of the lower indices. The general trend is however clear. One gets a relation 
of the form 

( fjfj+nfj+m )' = Y ( fj-afj+n-a' fj+m-a" ) + correction terms. (68) 



z 3 

a,a'a" 



The correction terms refer to the cases where the indices j, j + n and j + m are either 
equal or nearest neighbors. 

Consider now the class of uncorrelated distributions which coincide with the sta- 



tionary state flTT|) up to the second moment but start to deviate at the level of the third 
moment. The difference of the third order correlations with respect to the stationary 
state then only has a value for n = m = 0. Such a disturbance gradually spreads over 
the 2(d — 1) dimensional space of the indices n and m. The decay rate is therefore 
as l^ 1 '^ which is faster than the /( 1-ci )/ 2 of the slowest mode. The correction terms 
in ( j68|) mildly modify the spreading of the correlations described by the first term on 
the right hand side of (|68"D . They do no change the power of the decay but change the 
amplitude in front of the power. 

The same story holds for the construction of the higher modes. By further differ- 
entiation of relation Q58| ) one derives the transformation law of the functions which 
have 3 factors of the type Si/(z+ Sj) extra over the stationary state. Such functions do 
not modify the normalisation, the mean force and the second moment and therefore 
couple to the decay of the third moment as we have described above. On the level of 
the single-site distribution the mode takes the form 



p(s) =p*{s) 



1 + ^0,0 



z + s 



(69) 



with the B§1 ~ 
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Figure 3: The slowest mode for the triangular packing at I = 50 and the fee packing 
at I = 10 



The generalization to higher moments and faster modes is evident. On the level of 
the single site distribution it simply means a change of the power of the last factor in 
(|69|). The coefficient in front of a power m decays as /( 1_d )( m_1 )/ 2 . 

The question that emerges is: can one observe these modes in the relaxation of an 
actual simulation? This question is related to the problem of whether an arbitrary 
initial distribution can be decomposed as a superposition of these modes. Already 
on the level of the single site distribution one observes that our modes do not form 
a complete set in the mathematical sense. All modes start with a power f^ z ~^ in 
the force distribution; hence the modes form a poor basis for the small / behavior. 
Nevertheless, the modes that we have constructed do describe the asymptotic decay. 
We should realize that before the asymptotic regime sets in, a fast process takes place. 
If we start out with a distribution with a finite probability density for small forces, 
this probability will be supressed rapidly in the region of small forces by the recursion 
relation, due to the available phase space. For a small resultant force /■ to occur, 
all the components qf_ a have to be small. Already after one recursion step, a finite 
probability density (in the single-site distribution) will be replaced by a probability 
density starting as f z ~ l (modulo logarithmic factors, which move to higher orders in the 
next recursion steps). Thus after a few steps, the decay of the probability distribution 
can be well described by the above constructed modes. 

In order to test the mode picture we have performed numerical simulations on a 
2-dimensional triangular lattice and a 3-dimensional fee lattice. We took the initial 
condition fi — 1 for each site i, which also means (ff) = 1. This initial distribution 
is a substantial deviation from the stationary state flIT|) . In Fig. || we have compared 
the observed single-site distributions with the analytical results. In the first figure, we 
show the single-site force distribution of the I = 50th layer in the triangular packing, 
after subtracting p*(f). The analytical curve is the slowest mode whose amplitude 



|(0) 



-1, 



has been computed from the recursive scheme (Q), with the initial value Aq 
corresponding to the initial second moment {ff ) — 1. As the comparison does not 
involve any free parameters the correspondence is remarkable; the more so since the 
deviations bear the signature of the next slowest mode (|69|). In the second figure, the 
simulation data is shown for the fee packing with the same initial distribution and 
/ = 10. In line with the fact that the modes decay faster in a higher dimension, the 
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agreement is even more impressive. 



8 Discussion 

We have calculated the decay towards the stationary state (|TTf) of initial force distribu- 
tions which have short-ranged force correlations and a sufficiently sharp distribution of 



the total force (19). Since we can fix the average force, the relaxation occurs for the sec- 
ond and higher order force correlations. The short-ranged correlations become longer 
ranged while their amplitudes diminish in a diffusion process in correlation space. This 
accounts for the slow algebraic relaxation of the second moments of the force. These 
moments also dictate the slowest modes of relaxation of the full distribution function. 
For arbitrary initial conditions we find the asymptotic form fl66[) for the relaxation 
of the single-site distribution. We have focussed on the triangular and the fee pack- 
ing as examples, but the results can be easily extended to other packings and other 
connections between layers. The number of connections z determines the form of the 
stationary state and the decay modes. The dimension of the system detemines the 
relaxation powers of the modes. 

The stationary state ( |TTD is not the only stationary state as can already be seen 
from the fact that each initial distribution of the total force stays invariant. If we 
exclude, however, long-ranged correlations in the initial state and restrict ourselves to 
sharp distributions of the total force in the sense of (fl9|) , the stationary state can only 
marginally differ from ( |TT| ) as is explained in the appendix. 

It is not remarkable that the product stationary state ([ll]) results also from a 
mean-field approximation. However, it is interesting that the single-site distributions 
(|66| ) etc. are also found as modes of the mean-field approximation. Yet it would be 
misleading to conclude that the mean-field approximation is accurate. The mean-field 
modes decay exponentially as a consequence of ignoring correlations in the system. An 
initially uncorrelated state develops correlations which get longer in range and smaller 
in amplitude. The long persistence of these correlations explains why the mean-field 
approximation does not accurately describe the relaxation. 

In our calculations we have relied heavily on the mathematical simplifications due 
to the uniform distribution of the q variables. The slow approach to the stationary 
state is a robust property which holds for arbitrary q distributions. In fact it is easy to 
show that the decay of the second moments for an arbitrary q distribution is similar to 
what has been found in Section [| (see also [[J). However, the identity Q55| ) only holds 
for the uniform distribution. So there is as yet no general procedure to construct the 
explicit form of the stationary state and the slow modes. Coppersmith et al. || have 
found a countable set of q distributions for which the stationary state is uncorrelated, 
on the basis of a generalization of (0). For these distributions the construction of the 
modes is completely similar to what we have done for the uniform distribution. In 
general, the asymptotic distribution will exhibit correlations, albeit in higher moments 
than the second. The question how the correlations in the asymptotic state are related 
to the properties of the q distribution remains intriguing, even if the correlations are 
of higher order and likely to be small in magnitude. 
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A Other Stationary States 

There are two ways to investigate the occurrence of a stationary state. We can do this 
on the level of the moments or of the full probability distribution. As the moments are 
simpler we consider first the second moment. We take the total force distribution to 
be sharp according to fll~9|). Imposing this condition on the stationary state (|33| ) gives 
for the coefficient C 

So indeed the stationary values ( |33| ) approach those of ([11]) in the limit N — > oo. 

Similar arguments can be given for the stationary state fl57j) following from the scale 
invariance of the transformation. This state does not correspond to a good distribution 
as its norm is zero. This problem can be avoided by adding P*(s) and allowing for a 
different scale, i.e. 

*® = ( 1+A ?Trk)- (71) 

Unlike the slowest mode, this particular addition to P*(s) does affect the average force. 
Rescaling the forces such that (fi) = l, requires \z = 1 — A. The second moments of 
d?|) are 

(f l f l+n ) = (l-A*) Z -±^. (72) 

z 

Thus the distribution is sharp in the sense of Eq. fll9]) if 



A 2 ~ (l/z - c)/N, (73) 

which means that A is 0(1/ \fN) (or for c = l/z corresponding to the stationary 
state (ITT)) . Thus the stationary state (fn]) becomes equal to the stationary state (|TTD 



in the thermodynamic limit. The long-ranged correlations are indeed of the order 1/N. 

The same conclusion can be drawn with respect to the stationary state correspond- 
ing to Q8"5D which reads 

P*(s) = P*(s) (l + B J2Qn(s) + §Qo(*)) • (74) 
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The allowed values for B follow again from computing (F 2 ) for the distribution ( |74|) 
and equating this with the expression (|19D 

(l + 2B/z 2 )(N 2 + N/z) = N 2 + cN or B ~ z 2 (c - l/z)/2N. (75) 

So B turns out to be 0(1/ N) showing that the stationary state implied by ( |6~5"D is 
indeed only a marginal deviation from the stationary state (|TT|). It has long-ranged 
correlations of the order of 1/N. Similar arguments can be given for the other stationary 
states following from the scale invariance of the transformation. For finite systems such 
states are different, but in the thermodynamic limit they all merge in the stationary 
state QTTD if we insist on the sharpness of the total force distribution. 
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